Noncollapsibility of the odds ratio

code
analysis
IPW-analysis
Author

André Meichtry

Published

December 2, 2023

Noncollapsibility

Assume having a binary outcome \(Y\) and a treatment indicator \(X\) denoting if patient was randomised to treatment (\(X=1\)) or control (\(X=0\)). The analysis would be to fit a simple logistic regression model

\[ P(Y=1|X)=\text{expit}(\alpha_0+\alpha_1X), \]

with \(\text{expit}(\alpha_0+\alpha_1X)=\frac{\exp(\alpha_0+\alpha_1X)}{1+\exp(\alpha_0+\alpha_1X)}\).

Often, we write the model with the linear predictor, that is, the logit transformed expectation,

\[ \boxed{\text{logit}P(Y=1|X)=\alpha_0+{\color{red}\alpha_1}X}. \]

Of course, there will be other factors that influence the probability that \(Y=1\). Assume another variable \(C\) which influences the probability that \(Y=1\). The conditional model would be

\[ \boxed{\text{logit}P(Y=1|X,C)=\beta_0+{\color{red}\beta_1}X+\beta_2C}. \]

Now, under randomisation, \(C\) and \(X\) are (population) independent – that is – there is no confounding of the effect of \(X\) on \(Y\) by \(C\) in place.

Though \(C\) is not a confounder, in general, \({\color{red}\alpha_1 \neq \beta_1}\). The parameters (log odds ratios) in the marginal and conditional model are different. The proof is based on William’s Tower rule for conditional expectations:

\[ \begin{align} E(Y|X)=&E(E(Y|X,C)|X)\\ =&E(\text{expit}(\beta_0+\beta_1X+\beta_2C)|X)\\ \neq &\text{expit}(\beta_0+\beta_1X+\beta_2E(C|X)) \end{align} \]

With other link functions (in log-linear models) and identity (in linear models)), however, the two quantities are equal.

Example

We simulate data from logistic GLM with non-confounding covariate \(C\):

set.seed(3)
N <- 400
C <- sort(runif(N,-8,8)) #non-confounding C
X<-sample(c(0,1),N,replace=TRUE)
beta0<-0
beta1<-log(10)
beta2<-1
##Logistic 
etai<-beta0+beta1*X+beta2*C
pii<-exp(etai)/(1+exp(etai))
Ydich <-rbinom(X,size=1,prob=pii)
dat<-data.frame(C,X,Ydich)
psych::headTail(dat)
C X Ydich
1 -7.87 0 0
2 -7.87 0 0
3 -7.86 0 0
4 -7.81 0 0
397 7.93 1 1
398 7.94 0 1
399 7.95 0 1
400 7.95 0 1
modC<-glm(Ydich~X+C,data=dat,family="binomial")
modM<-glm(Ydich~X,data=dat,family="binomial")

Marginal model

equatiomatic::extract_eq(modM)

\[ \log\left[ \frac { P( \operatorname{Ydich} = \operatorname{1} ) }{ 1 - P( \operatorname{Ydich} = \operatorname{1} ) } \right] = \alpha + \beta_{1}(\operatorname{X}) \]

Conditional model

equatiomatic::extract_eq(modC)

\[ \log\left[ \frac { P( \operatorname{Ydich} = \operatorname{1} ) }{ 1 - P( \operatorname{Ydich} = \operatorname{1} ) } \right] = \alpha + \beta_{1}(\operatorname{X}) + \beta_{2}(\operatorname{C}) \]

Illustrate non-collapsibility

pred<-predict(modM,type="response",se.fit=TRUE)
predgA<-predict(modC,newdata=data.frame(C=C[X==0],X=0),se.fit=TRUE,type="response")
predgB<-predict(modC,newdata=data.frame(C=C[X==1],X=1),se.fit=TRUE,type="response")
plot(C,Ydich,col=c("blue","red")[X+1],ylab="Probability")
lines(sort(C[X==0]),sort(predgA$fit),col="blue")
lines(sort(C[X==1]),sort(predgB$fit),col="red")
abline(a=mean(pii[X==0]),b=0,col="blue",lty=2)
abline(a=mean(pii[X==1]),b=0,col="red",lty=2)
Log odds ratio for group is larger when we adjust for \(C\) (Differences in logits between solid lines relative to distance between dashed lines)

Marginal model

summary(modM)$coef
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.08004     0.1415 -0.5655 0.571710
X            0.63377     0.2040  3.1071 0.001889

Conditional model

Effect of \(X\) is larger, even though \(C\) is not a confounder. This is due to non-collapsibility of the odds ratio

summary(modC)$coef
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) -0.05395     0.3208 -0.1682 8.665e-01
X            2.46270     0.5209  4.7280 2.267e-06
C            1.04913     0.1200  8.7453 2.223e-18

IPW-model

Estimate inverse probability weights to fit marginal structural models in a point treatment situation.

Weights

library(ipw)
temp <- ipwpoint(exposure=X,family = "binomial",link = "logit",data=dat,numerator =~1,denominator = ~ C)
summary(temp$ipw.weights)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.989   0.994   1.000   1.000   1.006   1.012 
dat$sw<-temp$ipw.weights

Marginal structural model

Marginal structural model for the causal effect of \(X\) on \(Y\) corrected for confounding by \(C\) using inverse probability weighting

modW<-glm(Ydich~X,weights=sw,data=dat,family=quasibinomial)
summary(modW)$coef
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.0683     0.1419 -0.4814 0.630474
X             0.6117     0.2043  2.9938 0.002927
Loading required package: survey
Loading required package: grid
Loading required package: Matrix
Loading required package: survival

Attaching package: 'survey'
The following object is masked from 'package:graphics':

    dotchart
msm <- (svyglm(Ydich ~ X, family="binomial",design = svydesign(~ 1, weights =~sw,data = dat)))
Warning in eval(family$initialize): non-integer #successes in a binomial glm!
summary(msm)

Call:
svyglm(formula = Ydich ~ X, design = svydesign(~1, weights = ~sw, 
    data = dat), family = "binomial")

Survey design:
svydesign(~1, weights = ~sw, data = dat)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.0683     0.1417   -0.48   0.6301
X             0.6117     0.2042    3.00   0.0029

(Dispersion parameter for binomial family taken to be 1.003)

Number of Fisher Scoring iterations: 4

Comparison

modelsummary::modelsummary(models=list("marginal"=modM,"conditional"=modC,"IPW"=modW))
tinytable_go2gskbe8m18o691eacw
marginal conditional IPW
(Intercept) -0.080 -0.054 -0.068
(0.142) (0.321) (0.142)
X 0.634 2.463 0.612
(0.204) (0.521) (0.204)
C 1.049
(0.120)
Num.Obs. 400 400 400
AIC 543.4 152.6
BIC 551.4 164.5
Log.Lik. -269.718 -73.275
F 9.654 38.281 8.963
RMSE 0.49 0.23 0.49
Important
  • conditional-marginal=noncollaps+confouding
  • IPW-marginal=confounding
  • conditional-IPW=noncollaps

Smoking and children weight

str(bwt)
'data.frame':   189 obs. of  9 variables:
 $ low     : int  0 0 0 0 0 0 0 0 0 0 ...
 $ age     : int  19 33 20 21 18 21 22 17 29 26 ...
 $ lwt     : int  182 155 105 108 107 124 118 103 123 113 ...
 $ smoke   : int  0 0 1 1 1 0 0 0 1 1 ...
 $ ht      : int  0 0 0 0 0 0 0 0 0 0 ...
 $ ui      : int  1 0 0 1 1 0 0 0 0 0 ...
 $ ftv     : int  0 3 1 2 0 0 1 1 1 0 ...
 $ bwt     : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
 $ socclass: Factor w/ 3 levels "I","II","III": 2 3 1 1 1 3 1 3 1 1 ...
glm_smoke <- glm(low ~ smoke, data = bwt, family = binomial)
summary(glm_smoke)$coef
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  -1.0871     0.2147  -5.062 4.142e-07
smoke         0.7041     0.3196   2.203 2.762e-02
glm_multi_1 <- glm(low ~ smoke + socclass, data = bwt, family = binomial)
summary(glm_multi_1)$coef
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)   -1.841     0.3529  -5.216 1.828e-07
smoke          1.116     0.3692   3.023 2.507e-03
socclassII     1.084     0.4900   2.212 2.693e-02
socclassIII    1.109     0.4003   2.769 5.618e-03
temp <- ipwpoint(exposure=smoke,family = "binomial",link = "logit",data=bwt,
                 numerator=~1,denominator=~socclass)
summary(temp$ipw.weights)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.723   0.723   0.741   1.000   1.328   2.186 
bwt$sw<-temp$ipw.weights
modW<-glm(low~smoke,weights=sw,data=bwt,family=quasibinomial)
summary(modW)$coef
            Estimate Std. Error t value  Pr(>|t|)
(Intercept)   -1.277     0.2270  -5.625 6.696e-08
smoke          0.937     0.3283   2.854 4.798e-03
Important
  • conditional-marginal=noncollaps+confouding
  • IPW-marginal=confounding
  • conditional-IPW=noncollaps
modelsummary::modelsummary(models=list("marginal"=glm_smoke,"conditional"=glm_multi_1,"IPW"=modW))
tinytable_mok38oij1jh1rmtbb659
marginal conditional IPW
(Intercept) -1.087 -1.841 -1.277
(0.215) (0.353) (0.227)
smoke 0.704 1.116 0.937
(0.320) (0.369) (0.328)
socclassII 1.084
(0.490)
socclassIII 1.109
(0.400)
Num.Obs. 189 189 189
AIC 233.8 228.0
BIC 240.3 240.9
Log.Lik. -114.902 -109.987
F 4.852 4.294 8.148
RMSE 0.46 0.45 0.46

The view of potential outcomes

Simulate some data

library(simstudy)
# define the data
defB <- defData(varname = "C", formula =0.27, 
                dist = "binary")
defB <- defData(defB, varname = "Y0", formula = "-2.5 + 1.75*C", 
                dist = "binary", link = "logit")
defB <- defData(defB, varname = "Y1", formula = "-1.5 + 1.75*C", 
                dist = "binary", link = "logit")
defB <- defData(defB, varname = "X", formula = "0.315 + 0.352 *C", 
                dist = "binary")
defB <- defData(defB, varname = "Y", formula = "Y0 + X * (Y1 - Y0)", 
                dist = "nonrandom")
defB
varname formula variance dist link
C 0.27 0 binary identity
Y0 -2.5 + 1.75*C 0 binary logit
Y1 -1.5 + 1.75*C 0 binary logit
X 0.315 + 0.352 *C 0 binary identity
Y Y0 + X * (Y1 - Y0) 0 nonrandom identity
set.seed(2002)
Nnew<-1000
dtB <- genData(Nnew, defB)
dtB
id C Y0 Y1 X Y
1 0 0 0 0 0
2 0 0 1 0 0
3 1 0 0 1 0
4 0 0 0 0 0
5 1 1 1 0 1
6 1 0 0 1 0
7 1 0 0 1 0
8 1 0 1 0 0
9 1 1 1 1 1
10 1 0 0 1 0
11 0 0 1 1 1
12 0 0 0 1 0
13 1 0 0 0 0
14 0 0 1 0 0
15 0 0 1 1 1
16 0 0 0 1 0
17 0 0 0 0 0
18 0 0 1 0 0
19 0 0 0 0 0
20 0 0 0 0 0
21 0 0 0 0 0
22 1 1 0 1 0
23 0 0 0 0 0
24 0 0 1 0 0
25 0 0 0 0 0
26 0 0 0 0 0
27 0 0 0 1 0
28 0 1 1 0 1
29 0 0 0 0 0
30 0 0 0 0 0
31 1 0 1 1 1
32 0 0 0 0 0
33 1 1 1 0 1
34 0 0 0 0 0
35 0 0 1 0 0
36 0 0 0 1 0
37 0 0 0 1 0
38 0 0 0 0 0
39 0 0 0 1 0
40 1 1 0 1 0
41 1 0 0 0 0
42 1 1 0 1 0
43 0 0 0 0 0
44 1 1 1 1 1
45 0 0 0 1 0
46 0 0 0 1 0
47 0 0 0 1 0
48 0 0 0 0 0
49 0 0 1 0 0
50 0 0 0 0 0
51 0 0 0 0 0
52 0 0 0 1 0
53 0 0 0 0 0
54 0 0 0 1 0
55 0 0 0 1 0
56 0 0 0 0 0
57 0 0 0 0 0
58 1 0 1 0 0
59 0 0 0 1 0
60 1 0 1 1 1
61 0 0 0 0 0
62 1 1 1 1 1
63 0 0 1 1 1
64 0 0 0 0 0
65 1 0 1 1 1
66 1 0 0 1 0
67 0 0 1 0 0
68 1 1 0 0 1
69 0 0 0 0 0
70 0 0 0 0 0
71 0 0 0 0 0
72 0 0 0 0 0
73 0 0 0 0 0
74 1 1 1 1 1
75 0 0 0 0 0
76 0 0 0 0 0
77 0 0 0 0 0
78 0 0 1 1 1
79 0 0 1 0 0
80 0 0 0 0 0
81 0 0 1 0 0
82 0 0 0 0 0
83 0 0 0 0 0
84 1 0 1 1 1
85 0 0 1 0 0
86 1 0 1 1 1
87 0 0 1 0 0
88 0 0 1 0 0
89 0 0 1 0 0
90 0 0 0 0 0
91 0 0 0 0 0
92 1 0 0 1 0
93 1 0 1 0 0
94 0 0 0 0 0
95 1 0 1 0 0
96 0 1 0 1 0
97 1 0 1 1 1
98 0 0 0 1 0
99 0 0 0 0 0
100 0 0 0 0 0
101 0 0 0 0 0
102 0 0 0 1 0
103 1 0 0 1 0
104 1 1 1 1 1
105 0 0 0 0 0
106 0 0 0 0 0
107 0 0 0 0 0
108 0 1 0 1 0
109 0 0 0 0 0
110 0 1 0 0 1
111 0 0 1 0 0
112 1 0 0 1 0
113 0 0 0 1 0
114 0 0 0 0 0
115 0 0 0 0 0
116 0 0 0 0 0
117 0 0 0 1 0
118 0 0 1 0 0
119 0 0 0 1 0
120 0 0 0 0 0
121 1 0 1 1 1
122 1 0 0 1 0
123 0 0 0 0 0
124 0 0 0 0 0
125 0 0 1 0 0
126 1 0 1 0 0
127 1 1 0 1 0
128 0 0 0 0 0
129 1 0 1 1 1
130 0 0 0 1 0
131 0 0 1 1 1
132 1 0 0 1 0
133 1 1 0 1 0
134 0 0 0 0 0
135 0 0 1 1 1
136 0 0 0 0 0
137 0 0 0 0 0
138 0 0 0 0 0
139 0 1 0 0 1
140 0 0 0 0 0
141 0 0 0 1 0
142 1 0 1 1 1
143 0 0 0 0 0
144 1 0 0 0 0
145 0 0 0 1 0
146 0 0 0 1 0
147 0 0 0 1 0
148 0 0 0 0 0
149 0 0 0 0 0
150 0 0 0 0 0
151 0 0 0 0 0
152 0 0 1 0 0
153 1 0 1 1 1
154 0 0 0 0 0
155 0 0 1 0 0
156 1 1 0 1 0
157 0 1 0 0 1
158 0 0 0 0 0
159 0 0 0 0 0
160 0 0 1 0 0
161 0 0 1 0 0
162 1 0 0 0 0
163 0 0 1 0 0
164 1 0 0 1 0
165 0 0 0 0 0
166 0 0 0 0 0
167 0 0 0 0 0
168 0 0 0 0 0
169 1 1 0 1 0
170 0 0 1 1 1
171 0 0 0 0 0
172 0 0 0 1 0
173 1 1 0 1 0
174 1 0 1 1 1
175 0 0 0 0 0
176 0 0 0 0 0
177 0 0 0 1 0
178 1 0 0 1 0
179 0 0 0 1 0
180 1 1 0 1 0
181 1 1 1 0 1
182 0 0 0 1 0
183 0 0 0 1 0
184 1 0 1 0 0
185 0 0 0 1 0
186 0 0 0 0 0
187 0 0 0 1 0
188 0 0 0 1 0
189 0 0 0 0 0
190 0 0 0 0 0
191 0 0 0 0 0
192 0 0 0 0 0
193 0 0 0 0 0
194 0 0 0 0 0
195 0 0 0 1 0
196 1 0 1 1 1
197 0 0 0 0 0
198 1 0 1 1 1
199 1 0 0 1 0
200 0 0 0 1 0
201 0 1 0 1 0
202 1 0 0 0 0
203 1 1 1 0 1
204 0 0 0 0 0
205 0 0 0 0 0
206 0 0 0 0 0
207 0 0 1 0 0
208 1 0 1 0 0
209 0 0 1 1 1
210 1 0 0 1 0
211 0 0 0 0 0
212 0 0 0 1 0
213 0 0 0 0 0
214 1 1 0 1 0
215 0 0 1 1 1
216 0 0 0 0 0
217 1 0 0 1 0
218 0 0 0 1 0
219 0 0 0 0 0
220 0 0 0 1 0
221 0 0 0 0 0
222 0 0 0 0 0
223 1 1 0 1 0
224 0 0 0 0 0
225 1 0 1 0 0
226 0 0 0 0 0
227 0 0 0 1 0
228 1 0 0 0 0
229 0 0 0 0 0
230 1 0 1 1 1
231 0 0 0 1 0
232 0 0 0 0 0
233 0 0 0 0 0
234 1 0 1 1 1
235 0 0 1 0 0
236 1 0 1 0 0
237 0 0 0 0 0
238 1 1 0 1 0
239 0 0 0 0 0
240 0 0 0 0 0
241 0 0 0 0 0
242 0 0 1 0 0
243 0 0 0 0 0
244 0 0 0 0 0
245 0 0 0 0 0
246 1 0 1 1 1
247 1 0 1 0 0
248 1 0 0 0 0
249 0 0 0 0 0
250 0 0 0 0 0
251 1 0 0 1 0
252 0 0 0 1 0
253 0 0 0 0 0
254 0 0 0 0 0
255 0 0 0 1 0
256 0 0 0 0 0
257 0 0 0 0 0
258 1 1 1 0 1
259 0 0 0 0 0
260 0 0 0 0 0
261 0 0 0 0 0
262 0 0 0 0 0
263 0 0 0 1 0
264 0 0 1 1 1
265 0 0 0 0 0
266 1 0 0 1 0
267 0 0 0 0 0
268 0 0 0 0 0
269 0 0 0 0 0
270 0 0 0 0 0
271 0 0 1 1 1
272 0 0 0 0 0
273 0 0 0 0 0
274 1 0 1 0 0
275 1 0 1 0 0
276 0 0 0 0 0
277 0 0 0 1 0
278 0 0 0 0 0
279 0 0 0 1 0
280 0 0 0 1 0
281 1 1 1 0 1
282 0 0 0 0 0
283 0 0 0 0 0
284 1 0 1 0 0
285 1 1 1 0 1
286 0 0 0 0 0
287 0 0 0 0 0
288 1 0 1 0 0
289 0 0 0 1 0
290 0 0 0 0 0
291 0 0 1 0 0
292 0 0 0 1 0
293 1 0 1 1 1
294 0 0 0 0 0
295 1 0 1 1 1
296 0 0 0 1 0
297 0 0 0 1 0
298 1 1 1 1 1
299 0 0 0 1 0
300 0 0 1 0 0
301 0 0 1 1 1
302 1 0 1 1 1
303 0 0 0 0 0
304 0 0 0 0 0
305 1 0 0 1 0
306 0 0 1 0 0
307 0 0 1 0 0
308 1 0 0 1 0
309 1 0 1 1 1
310 0 0 0 1 0
311 0 1 0 0 1
312 0 0 1 0 0
313 0 0 0 0 0
314 0 0 0 1 0
315 0 0 0 1 0
316 0 0 0 0 0
317 1 0 1 0 0
318 0 1 0 0 1
319 1 0 1 1 1
320 1 1 1 1 1
321 0 0 0 0 0
322 0 0 1 1 1
323 0 0 0 0 0
324 0 0 0 0 0
325 1 0 1 1 1
326 0 0 0 0 0
327 0 0 0 1 0
328 0 0 0 1 0
329 0 0 0 0 0
330 0 0 0 0 0
331 0 0 0 0 0
332 0 0 1 0 0
333 0 0 0 0 0
334 0 0 0 0 0
335 0 0 0 1 0
336 0 0 0 1 0
337 0 0 0 0 0
338 0 0 0 0 0
339 0 0 0 1 0
340 0 0 0 0 0
341 0 0 1 0 0
342 1 0 1 0 0
343 0 0 0 0 0
344 0 0 0 0 0
345 1 0 0 1 0
346 1 0 0 1 0
347 1 1 1 1 1
348 0 0 0 1 0
349 0 0 0 0 0
350 0 0 0 0 0
351 0 0 0 1 0
352 0 0 0 0 0
353 1 0 0 1 0
354 0 0 0 0 0
355 1 1 0 1 0
356 0 0 0 0 0
357 1 1 0 1 0
358 1 0 1 1 1
359 1 0 0 1 0
360 0 0 1 0 0
361 0 0 0 0 0
362 1 0 1 0 0
363 0 0 0 0 0
364 0 0 0 1 0
365 1 0 1 1 1
366 0 0 0 0 0
367 0 0 0 0 0
368 1 0 1 1 1
369 0 0 0 0 0
370 0 0 0 0 0
371 0 1 0 1 0
372 0 0 0 1 0
373 0 0 0 1 0
374 0 0 0 0 0
375 1 0 1 0 0
376 0 0 1 1 1
377 1 1 0 1 0
378 0 0 0 0 0
379 0 0 0 1 0
380 1 0 0 1 0
381 1 0 1 1 1
382 1 0 1 1 1
383 0 0 0 0 0
384 0 0 1 0 0
385 0 0 0 0 0
386 0 1 0 1 0
387 1 0 0 0 0
388 0 0 0 1 0
389 0 0 0 0 0
390 0 0 0 0 0
391 1 0 1 1 1
392 1 0 0 1 0
393 0 0 0 0 0
394 0 0 0 1 0
395 0 0 0 0 0
396 1 0 0 1 0
397 0 0 0 1 0
398 0 0 0 0 0
399 0 0 0 0 0
400 0 0 0 0 0
401 1 1 1 1 1
402 1 1 1 0 1
403 1 0 1 1 1
404 1 0 1 1 1
405 0 0 0 1 0
406 0 0 0 0 0
407 1 0 1 1 1
408 0 0 0 1 0
409 0 0 1 0 0
410 1 1 1 1 1
411 0 0 0 0 0
412 1 0 1 1 1
413 0 0 0 0 0
414 0 0 0 0 0
415 0 1 1 0 1
416 1 0 1 0 0
417 0 0 0 1 0
418 0 0 0 1 0
419 0 0 1 0 0
420 0 0 1 1 1
421 1 0 0 1 0
422 0 0 0 1 0
423 0 0 0 0 0
424 0 0 0 0 0
425 0 0 0 1 0
426 0 0 0 1 0
427 0 0 0 1 0
428 0 0 1 0 0
429 0 0 0 0 0
430 1 0 1 1 1
431 0 0 0 0 0
432 0 0 0 0 0
433 0 0 1 1 1
434 0 0 0 1 0
435 0 0 0 1 0
436 0 0 0 1 0
437 0 0 0 1 0
438 0 0 1 1 1
439 0 0 0 0 0
440 0 0 0 0 0
441 1 0 0 1 0
442 0 0 0 0 0
443 1 0 0 0 0
444 0 0 0 0 0
445 0 0 0 0 0
446 0 0 1 1 1
447 1 1 0 0 1
448 0 0 0 0 0
449 0 0 0 0 0
450 0 0 0 0 0
451 1 0 0 1 0
452 1 0 0 0 0
453 0 0 0 0 0
454 1 0 0 1 0
455 0 0 0 1 0
456 0 0 0 0 0
457 0 0 0 0 0
458 0 0 0 0 0
459 0 0 0 0 0
460 0 0 0 1 0
461 1 0 1 1 1
462 1 0 1 1 1
463 0 0 0 0 0
464 0 0 0 1 0
465 0 0 0 0 0
466 0 0 0 0 0
467 1 1 1 1 1
468 0 1 0 0 1
469 1 1 0 1 0
470 0 0 0 1 0
471 0 0 1 0 0
472 0 0 0 1 0
473 0 0 0 0 0
474 1 1 0 1 0
475 0 0 1 0 0
476 0 0 0 0 0
477 0 1 0 0 1
478 0 0 0 1 0
479 1 0 0 0 0
480 0 0 0 0 0
481 0 0 1 1 1
482 0 0 0 1 0
483 1 0 0 1 0
484 0 0 0 1 0
485 0 0 0 0 0
486 0 0 0 0 0
487 0 0 0 1 0
488 1 0 1 0 0
489 0 0 0 1 0
490 0 0 0 0 0
491 0 0 0 0 0
492 1 0 0 1 0
493 0 0 0 1 0
494 0 0 0 1 0
495 1 1 1 1 1
496 0 0 0 0 0
497 1 1 1 1 1
498 0 0 0 0 0
499 0 0 0 1 0
500 1 0 1 1 1
501 0 0 0 1 0
502 0 1 0 1 0
503 0 0 0 1 0
504 0 0 0 0 0
505 1 0 0 1 0
506 1 1 1 1 1
507 0 0 0 0 0
508 0 0 0 0 0
509 0 0 1 1 1
510 1 0 1 0 0
511 1 1 1 1 1
512 0 0 0 0 0
513 0 0 0 0 0
514 0 0 0 0 0
515 1 0 0 0 0
516 0 0 0 1 0
517 0 0 0 1 0
518 1 0 0 0 0
519 1 1 1 1 1
520 0 0 0 0 0
521 0 0 0 0 0
522 0 0 0 1 0
523 0 0 0 0 0
524 0 0 0 1 0
525 1 1 0 1 0
526 0 0 0 1 0
527 0 0 1 1 1
528 0 0 0 1 0
529 0 0 1 0 0
530 0 0 0 0 0
531 0 0 0 0 0
532 0 0 0 1 0
533 1 0 1 1 1
534 0 0 1 1 1
535 0 0 0 1 0
536 0 0 0 1 0
537 0 0 0 0 0
538 0 0 0 0 0
539 0 0 0 0 0
540 0 0 0 1 0
541 0 0 0 0 0
542 0 0 0 1 0
543 1 1 1 1 1
544 1 1 0 1 0
545 0 0 1 0 0
546 0 0 0 1 0
547 0 0 0 1 0
548 0 0 0 1 0
549 0 0 0 0 0
550 0 0 0 0 0
551 0 0 0 0 0
552 0 0 0 1 0
553 1 0 1 1 1
554 1 0 1 1 1
555 1 0 1 0 0
556 0 0 0 0 0
557 0 0 0 0 0
558 0 0 0 0 0
559 1 0 1 0 0
560 1 0 1 0 0
561 1 0 0 1 0
562 0 0 0 0 0
563 0 0 0 0 0
564 0 0 0 0 0
565 1 0 1 1 1
566 0 0 0 1 0
567 0 0 1 1 1
568 0 0 0 0 0
569 0 0 0 0 0
570 0 0 1 0 0
571 0 0 0 1 0
572 0 0 1 0 0
573 0 0 0 1 0
574 1 0 1 0 0
575 0 0 0 1 0
576 0 0 0 1 0
577 1 0 0 1 0
578 0 0 0 0 0
579 0 0 1 0 0
580 1 0 0 0 0
581 0 0 1 0 0
582 0 0 1 1 1
583 0 0 1 1 1
584 0 0 0 0 0
585 1 0 1 0 0
586 0 0 0 0 0
587 0 0 1 0 0
588 0 0 0 0 0
589 0 1 0 1 0
590 0 0 0 0 0
591 0 0 0 1 0
592 0 0 0 0 0
593 0 1 0 1 0
594 1 0 1 1 1
595 1 0 1 0 0
596 0 1 1 1 1
597 1 0 1 1 1
598 0 0 0 0 0
599 1 0 1 0 0
600 1 0 1 1 1
601 1 0 0 1 0
602 0 0 0 0 0
603 0 0 1 1 1
604 1 1 0 0 1
605 0 0 0 0 0
606 0 0 0 1 0
607 0 0 0 1 0
608 0 1 0 0 1
609 1 0 1 0 0
610 0 0 0 0 0
611 0 0 0 1 0
612 0 0 1 0 0
613 0 0 0 1 0
614 0 0 1 0 0
615 0 0 0 0 0
616 0 0 0 1 0
617 0 0 0 1 0
618 0 0 0 1 0
619 0 0 0 0 0
620 1 0 0 1 0
621 0 0 0 0 0
622 0 0 0 0 0
623 0 0 1 0 0
624 0 0 1 1 1
625 0 0 0 0 0
626 1 1 0 1 0
627 0 1 0 1 0
628 1 1 1 1 1
629 1 1 0 1 0
630 0 0 0 1 0
631 1 0 1 1 1
632 0 0 0 1 0
633 0 0 0 1 0
634 0 0 0 0 0
635 0 0 0 0 0
636 0 0 0 0 0
637 0 0 0 1 0
638 0 1 0 0 1
639 1 0 0 1 0
640 0 0 0 0 0
641 0 0 0 0 0
642 0 1 0 0 1
643 0 0 0 1 0
644 0 0 0 1 0
645 0 0 0 0 0
646 1 1 0 1 0
647 0 0 0 0 0
648 0 0 1 0 0
649 1 1 1 1 1
650 0 0 0 1 0
651 0 0 0 0 0
652 0 0 0 0 0
653 1 0 1 0 0
654 0 0 0 0 0
655 0 0 0 0 0
656 0 0 0 0 0
657 0 0 0 1 0
658 0 0 0 0 0
659 1 0 1 0 0
660 0 0 1 0 0
661 0 0 0 1 0
662 0 0 0 1 0
663 0 0 0 0 0
664 1 1 1 0 1
665 0 1 0 1 0
666 0 0 0 1 0
667 0 0 0 1 0
668 1 0 0 1 0
669 0 0 0 0 0
670 0 0 0 0 0
671 0 1 0 0 1
672 0 0 0 0 0
673 0 1 0 0 1
674 0 0 0 1 0
675 0 0 0 0 0
676 1 0 0 0 0
677 0 0 0 0 0
678 0 1 0 0 1
679 0 1 0 1 0
680 0 0 0 0 0
681 0 0 0 1 0
682 1 0 1 1 1
683 0 1 0 0 1
684 0 0 0 0 0
685 0 0 0 0 0
686 0 0 0 0 0
687 0 0 0 1 0
688 0 0 0 0 0
689 0 0 0 1 0
690 1 0 1 0 0
691 0 0 0 1 0
692 0 0 0 0 0
693 1 0 0 1 0
694 0 0 1 0 0
695 0 0 0 0 0
696 0 0 0 1 0
697 0 0 1 0 0
698 0 0 0 0 0
699 0 0 0 0 0
700 0 0 0 0 0
701 0 0 1 0 0
702 0 0 0 1 0
703 1 1 1 1 1
704 0 0 1 1 1
705 0 0 0 0 0
706 0 0 0 1 0
707 0 0 0 1 0
708 0 0 0 0 0
709 0 0 0 0 0
710 0 0 0 0 0
711 0 0 1 0 0
712 0 0 1 0 0
713 1 0 1 0 0
714 0 0 0 0 0
715 0 0 0 0 0
716 0 0 1 0 0
717 0 0 0 0 0
718 0 0 1 0 0
719 1 0 0 1 0
720 0 1 1 0 1
721 0 0 0 0 0
722 0 0 1 0 0
723 0 0 0 1 0
724 0 0 0 1 0
725 0 0 0 0 0
726 0 0 0 1 0
727 0 0 0 0 0
728 0 0 0 0 0
729 0 0 0 0 0
730 1 0 0 1 0
731 0 0 0 0 0
732 0 0 0 0 0
733 0 0 0 0 0
734 0 0 1 0 0
735 0 0 0 0 0
736 0 0 0 0 0
737 0 0 0 0 0
738 0 0 0 1 0
739 0 0 1 0 0
740 0 0 0 0 0
741 1 0 0 0 0
742 0 0 0 0 0
743 0 0 0 0 0
744 0 0 1 0 0
745 0 0 1 1 1
746 0 0 0 0 0
747 0 0 0 1 0
748 0 0 1 1 1
749 0 1 0 1 0
750 0 0 0 0 0
751 0 0 0 1 0
752 0 0 0 0 0
753 0 0 1 1 1
754 0 0 0 0 0
755 0 0 0 0 0
756 0 0 1 0 0
757 1 1 1 1 1
758 0 0 1 0 0
759 1 0 1 1 1
760 0 0 0 1 0
761 0 1 0 0 1
762 1 0 1 0 0
763 0 1 0 0 1
764 0 0 0 1 0
765 0 0 0 0 0
766 0 0 0 0 0
767 0 0 0 0 0
768 1 0 1 1 1
769 0 0 0 1 0
770 0 0 0 1 0
771 1 0 1 1 1
772 0 0 0 1 0
773 0 0 0 0 0
774 1 1 1 0 1
775 0 0 0 1 0
776 1 1 0 1 0
777 0 0 0 0 0
778 0 0 1 1 1
779 0 0 0 0 0
780 1 1 0 0 1
781 0 1 0 0 1
782 1 1 0 1 0
783 0 0 0 0 0
784 0 0 0 0 0
785 1 0 1 1 1
786 1 1 1 0 1
787 0 1 0 0 1
788 0 0 0 1 0
789 1 0 0 1 0
790 0 0 0 1 0
791 0 0 0 1 0
792 1 1 0 0 1
793 0 0 0 0 0
794 0 0 0 0 0
795 0 1 0 1 0
796 0 0 1 0 0
797 0 0 0 0 0
798 1 0 0 1 0
799 1 1 1 0 1
800 0 0 0 0 0
801 0 0 0 1 0
802 0 0 1 0 0
803 0 0 0 1 0
804 0 0 0 0 0
805 0 0 0 0 0
806 0 0 0 0 0
807 0 0 0 0 0
808 0 0 0 0 0
809 1 1 0 1 0
810 0 0 0 0 0
811 0 1 0 1 0
812 0 0 0 0 0
813 1 0 0 0 0
814 1 1 1 1 1
815 0 0 0 0 0
816 1 0 1 1 1
817 0 1 1 1 1
818 1 0 0 1 0
819 1 0 0 1 0
820 0 0 0 0 0
821 0 0 1 0 0
822 0 0 0 0 0
823 0 1 1 0 1
824 0 0 0 0 0
825 0 0 0 1 0
826 0 0 0 0 0
827 0 0 1 1 1
828 0 1 0 0 1
829 0 0 0 0 0
830 1 1 1 1 1
831 1 0 0 1 0
832 1 1 1 0 1
833 0 0 0 0 0
834 0 0 0 0 0
835 0 0 0 0 0
836 0 0 0 0 0
837 0 0 0 0 0
838 0 0 1 0 0
839 0 0 1 1 1
840 0 0 0 0 0
841 0 1 1 1 1
842 0 0 0 1 0
843 1 1 1 1 1
844 0 0 0 1 0
845 1 0 1 1 1
846 0 0 0 0 0
847 1 0 1 0 0
848 0 0 0 0 0
849 1 0 1 0 0
850 0 0 0 0 0
851 0 0 0 0 0
852 0 0 0 0 0
853 1 0 0 1 0
854 1 0 0 0 0
855 0 0 0 0 0
856 0 0 0 0 0
857 0 0 0 1 0
858 1 1 0 0 1
859 0 0 0 1 0
860 0 0 0 0 0
861 0 0 0 1 0
862 0 1 0 0 1
863 1 0 1 1 1
864 0 0 0 1 0
865 0 0 0 0 0
866 0 0 1 0 0
867 0 1 0 0 1
868 1 0 1 1 1
869 1 0 1 0 0
870 0 0 0 0 0
871 1 0 1 0 0
872 0 0 1 0 0
873 0 0 0 0 0
874 0 0 0 0 0
875 0 0 0 0 0
876 0 0 0 0 0
877 1 0 0 0 0
878 0 1 1 0 1
879 0 0 0 0 0
880 1 0 1 0 0
881 0 0 0 0 0
882 0 0 0 0 0
883 0 0 0 0 0
884 0 0 0 1 0
885 1 0 1 1 1
886 0 0 0 1 0
887 1 0 0 1 0
888 1 0 1 1 1
889 0 0 0 0 0
890 1 1 0 0 1
891 1 0 1 1 1
892 0 0 1 0 0
893 1 0 1 0 0
894 0 0 0 0 0
895 0 0 0 0 0
896 1 1 0 1 0
897 0 0 0 0 0
898 0 0 0 0 0
899 0 0 0 0 0
900 0 0 0 1 0
901 1 0 0 1 0
902 1 0 0 1 0
903 1 1 0 1 0
904 1 1 0 0 1
905 0 0 1 1 1
906 1 0 1 1 1
907 0 0 0 0 0
908 0 0 1 0 0
909 0 0 1 0 0
910 0 0 0 0 0
911 0 0 1 1 1
912 1 1 1 1 1
913 0 0 1 0 0
914 1 0 1 1 1
915 0 0 0 1 0
916 1 0 1 0 0
917 0 0 0 0 0
918 0 0 0 1 0
919 0 0 0 1 0
920 0 0 0 0 0
921 0 0 0 0 0
922 0 0 0 0 0
923 0 0 0 0 0
924 0 0 0 1 0
925 0 0 0 1 0
926 0 0 1 0 0
927 0 0 0 0 0
928 1 0 0 0 0
929 0 0 0 0 0
930 0 0 0 0 0
931 1 0 0 1 0
932 0 0 0 0 0
933 1 0 1 1 1
934 0 0 0 1 0
935 0 0 1 1 1
936 1 1 1 1 1
937 0 0 0 0 0
938 0 1 1 1 1
939 0 0 0 1 0
940 0 0 0 1 0
941 0 0 0 0 0
942 0 0 0 0 0
943 1 0 1 0 0
944 1 1 1 1 1
945 0 1 0 1 0
946 0 0 1 1 1
947 1 0 1 1 1
948 0 0 0 0 0
949 0 0 0 0 0
950 1 0 1 1 1
951 1 1 0 1 0
952 0 0 0 0 0
953 1 0 1 0 0
954 1 1 0 0 1
955 0 0 1 1 1
956 0 0 0 0 0
957 0 0 0 0 0
958 0 0 1 1 1
959 0 0 0 0 0
960 0 1 0 1 0
961 0 0 0 0 0
962 0 0 0 0 0
963 0 0 1 1 1
964 0 0 1 0 0
965 0 0 1 0 0
966 0 0 0 0 0
967 0 0 1 1 1
968 1 0 0 1 0
969 0 0 0 0 0
970 0 0 0 1 0
971 0 0 0 0 0
972 0 0 0 0 0
973 0 0 0 0 0
974 0 1 0 0 1
975 0 0 0 1 0
976 0 0 0 1 0
977 1 0 0 1 0
978 0 0 0 0 0
979 1 0 1 1 1
980 1 0 0 1 0
981 0 0 0 1 0
982 0 0 0 1 0
983 0 0 0 0 0
984 0 0 0 0 0
985 0 0 0 0 0
986 1 1 1 0 1
987 0 0 0 0 0
988 0 0 0 0 0
989 0 0 0 1 0
990 0 0 0 0 0
991 1 1 0 1 0
992 0 0 0 0 0
993 0 1 0 0 1
994 1 1 0 1 0
995 0 0 0 1 0
996 0 0 1 0 0
997 0 0 0 1 0
998 1 0 1 1 1
999 0 0 0 0 0
1000 1 0 0 1 0

True causal effect (based on potential outcomes)

odds <- function (p) {
    return((p/(1 - p)))
}

dtB[, log( odds( mean(Y1) ) / odds( mean(Y0) ) )]
[1] 0.9494

Conditional effect

The true conditional causal effect of \(A\) is 1.

mc<-glm(Y ~ X + C , data = dtB, family="binomial")
mc

Call:  glm(formula = Y ~ X + C, family = "binomial", data = dtB)

Coefficients:
(Intercept)            X            C  
      -2.74         1.21         1.61  

Degrees of Freedom: 999 Total (i.e. Null);  997 Residual
Null Deviance:      964 
Residual Deviance: 797  AIC: 803

This estimate for \(X\) is a good estimate of the conditional effect in the population, based on the potential outcomes at each level of \(C\)

dtB[, .(LOR = log( odds( mean(Y1) ) / odds( mean(Y0) ) ) ), keyby = C]
C LOR
0 1.104
1 1.067

Marginal effect

The marginal estimate is biased both for the conditional effect and the marginal causal effect.

mm<-glm(Y ~ X , data = dtB, family="binomial")
mm

Call:  glm(formula = Y ~ X, family = "binomial", data = dtB)

Coefficients:
(Intercept)            X  
      -2.33         1.59  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      964 
Residual Deviance: 876  AIC: 880

Numerator and Denominator model

numModel <- glm(X ~ 1, data = dtB, family = "binomial")
denModel <- glm(X ~ C, data = dtB, family = "binomial")

Stabilized weights by hand

dtB[, pX0 := predict(numModel, type = "response")]
dtB[, pX := predict(denModel, type = "response")]
defB2 <- defDataAdd(varname = "IPW", 
                    formula = "(X*pX0+(1-X)*(1-pX0))/((X * pX) + ((1 - X) * (1 - pX)))", 
                    dist = "nonrandom")
dtB <- addColumns(defB2, dtB)
dtB[1:6]
id C Y0 Y1 X Y pX0 pX IPW
1 0 0 0 0 0 0.423 0.3347 0.8673
2 0 0 1 0 0 0.423 0.3347 0.8673
3 1 0 0 1 0 0.423 0.6718 0.6297
4 0 0 0 0 0 0.423 0.3347 0.8673
5 1 1 1 0 1 0.423 0.6718 1.7578
6 1 0 0 1 0 0.423 0.6718 0.6297
unique(dtB$IPW)
[1] 0.8673 0.6297 1.7578 1.2639

Stabilized weights with ipw package

tempsw<-ipw::ipwpoint(exposure=X,family="binomial",link="logit",numerator=~1,denominator = ~C,data=dtB)
tempw<-ipw::ipwpoint(exposure=X,family="binomial",link="logit",denominator = ~C,data=dtB)
unique(tempsw$ipw.weights)
[1] 0.8673 0.6297 1.7578 1.2639
mean(tempsw$ipw.weights)
[1] 1
unique(tempw$ipw.weights)
[1] 1.503 1.489 3.047 2.988
mean(tempw$ipw.weights)
[1] 2

Applying IPW

Unstabilized weights

summary(mw<-glm(Y ~ X , data = dtB, family="binomial", weights = tempw$ipw.weights))$coef
Warning in eval(family$initialize): non-integer #successes in a binomial glm!
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)   -2.069     0.1002 -20.645 1.076e-94
X              1.081     0.1229   8.801 1.354e-18

Stabilized weights

summary(mws<-glm(Y ~ X , data = dtB, family="binomial", weights = tempsw$ipw.weights))$coef
Warning in eval(family$initialize): non-integer #successes in a binomial glm!
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)   -2.069     0.1319 -15.682 1.996e-55
X              1.081     0.1713   6.312 2.760e-10

Comparison

modelsummary::modelsummary(models=list("marginal"=mm,"conditional"=mc,"IPW"=mw,"IPWs"=mws))
Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Warning in eval(family$initialize): non-integer #successes in a binomial glm!
tinytable_a1qfbesdj2xn5ksy0m40
marginal conditional IPW IPWs
(Intercept) -2.333 -2.736 -2.069 -2.069
(0.147) (0.165) (0.100) (0.132)
X 1.587 1.211 1.081 1.081
(0.180) (0.191) (0.123) (0.171)
C 1.614
(0.183)
Num.Obs. 1000 1000 1000 1000
AIC 880.1 802.5 1847.3 1004.4
BIC 889.9 817.3 1857.1 1014.2
Log.Lik. -438.049 -398.272 -921.632 -500.206
F 77.828 71.500 77.460 39.837
RMSE 0.37 0.35 0.37 0.37

We used R version 4.4.1 (R Core Team 2024) and the following R packages: ipw v. 1.2.1 (van der Wal and Geskus 2011), Matrix v. 1.7.0 (Bates, Maechler, and Jagan 2024), rio v. 1.0.1 (Chan et al. 2023), simstudy v. 0.7.1 (Goldfeld and Wujciak-Jens 2020), survey v. 4.4.2 (Lumley 2004, 2010, 2024), survival v. 3.6.4 [@].

Package Version Citation
base 4.4.1 R Core Team (2024)
ipw 1.2.1 van der Wal and Geskus (2011)
Matrix 1.7.0 Bates, Maechler, and Jagan (2024)
rio 1.0.1 Chan et al. (2023)
simstudy 0.7.1 Goldfeld and Wujciak-Jens (2020)
survey 4.4.2 Lumley (2004); Lumley (2010); Lumley (2024)
survival 3.6.4 @

Bibliography

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2024. Matrix: Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.
Chan, Chung-hong, Thomas J. Leeper, Jason Becker, and David Schoch. 2023. rio: A Swiss-Army Knife for Data File i/o. https://cran.r-project.org/package=rio.
Goldfeld, Keith, and Jacob Wujciak-Jens. 2020. simstudy: Illuminating Research Methods Through Data Generation.” Journal of Open Source Software 5 (54): 2763. https://doi.org/10.21105/joss.02763.
Lumley, Thomas. 2004. “Analysis of Complex Survey Samples.” Journal of Statistical Software 9 (1): 1–19.
———. 2010. Complex Surveys: A Guide to Analysis Using r: A Guide to Analysis Using r. John Wiley; Sons.
———. 2024. survey: Analysis of Complex Survey Samples.”
R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
van der Wal, Willem M., and Ronald B. Geskus. 2011. ipw: An R Package for Inverse Probability Weighting.” Journal of Statistical Software 43 (13): 1–23. https://doi.org/10.18637/jss.v043.i13.